library(tidyr)
library(dplyr)
library(plotly)
library(ggplot2)
library(lubridate)
library(stringr)
library(statsr)
library(ggmap)
library(graphics)

train <- read.csv("c:/Users/admin/Downloads/Pos/train/train.csv")

#Referências:
#https://www.novayork.net/taxi#
#http://carlosdelfino.eti.br/cursoarduino/geoprocessamento/calculando-distancias-com-base-em-coordenadas-de-gps/

Cálculando a distância Euclidiana em KM e em Milhas

d2r <- 0.017453292519943295769236

train <- mutate(train, DLA = (train$dropoff_latitude - train$pickup_latitude) * d2r,
                       DLO = (train$dropoff_longitude - train$pickup_longitude) * d2r)

train <- mutate(train, a = (sin(train$DLA / 2.0) ^ 2) + ((cos(train$pickup_latitude * d2r) ^ 2) * (sin(train$DLO / 2.0) ^ 2)))

train <- mutate(train, c = 2.0 * atan2(sqrt(train$a),sqrt(1.0 - train$a)))

train <- mutate(train, distEucKM = 6368.1 * train$c)

train$DLA <- NULL
train$DLO <- NULL
train$a <- NULL
train$c <- NULL

train <- mutate(train, distEucMilhas = 1.609 / train$distEucKM)

Calculando a distância de Manhattan em KM e em Milhas

train <- mutate(train, pontoTangLat = train$dropoff_latitude, pontoTangLong = train$pickup_longitude)

train <- mutate(train, DLAMan1 = (train$dropoff_latitude - train$pontoTangLat) * d2r,
                       DLOMan1 = (train$dropoff_longitude - train$pontoTangLong) * d2r)

train <- mutate(train, aMan1 = (sin(train$DLAMan1 / 2.0) ^ 2) + ((cos(train$pontoTangLat * d2r) ^ 2) * (sin(train$DLOMan1 / 2.0) ^ 2)))

train <- mutate(train, cMan1 = 2.0 * atan2(sqrt(train$aMan1),sqrt(1.0 - train$aMan1)))

train <- mutate(train, distManKM1 = 6368.1 * train$cMan1)

train$DLAMan1 <- NULL
train$DLOMan1 <- NULL
train$aMan1 <- NULL
train$cMan1 <- NULL

train <- mutate(train, DLAMan2 = (train$pickup_latitude - train$pontoTangLat) * d2r,
                       DLOMan2 = (train$pickup_longitude - train$pontoTangLong) * d2r)

train <- mutate(train, aMan2 = (sin(train$DLAMan2 / 2.0) ^ 2) + ((cos(train$pontoTangLat * d2r) ^ 2) * (sin(train$DLOMan2 / 2.0) ^ 2)))

train <- mutate(train, cMan2 = 2.0 * atan2(sqrt(train$aMan2),sqrt(1.0 - train$aMan2)))

train <- mutate(train, distManKM2 = 6368.1 * train$cMan2)

train$DLAMan2 <- NULL
train$DLOMan2 <- NULL
train$aMan2 <- NULL
train$cMan2 <- NULL

train <- mutate(train, distManKM = train$distManKM1 + train$distManKM2)
train <- mutate(train, distManMilhas = 1.609 / train$distManKM)

Calculando velocidade média

train <- mutate(train, Vmedia = train$distManKM / (train$trip_duration / 3600))

Calculando valor da corrida

train <- mutate(train, Pickuphms = str_sub(train$pickup_datetime, start = 11))

train <- mutate(train, valor = 2.5 + (1.56 * train$distManKM) + 0.5 + ifelse(train$Vmedia < 20,0.5 * (train$trip_duration / 60),0) + ifelse(hms(train$Pickuphms) >= hms('16:00:00') & hms(train$Pickuphms) <= hms('20:00:00'),1,0) + ifelse(hms(train$Pickuphms) >= hms('20:00:00') & hms(train$Pickuphms) <= hms('6:00:00'),0.5,0))

train <- mutate(train, valor = round(train$valor,2))

Configuração Inicial - Mapa de Calor

cfg <- list( x.min = -73.99525, x.max = -73.97545, x.step = 0.00018,
             y.min =  40.74545, y.max =  40.76525, y.step = 0.00018)

#Sequencias e quantidades de X e Y
cfg$x.lim <- seq( from = cfg$x.min, to = cfg$x.max, by = cfg$x.step)
cfg$y.lim <- seq( from = cfg$y.min, to = cfg$y.max, by = cfg$y.step)
cfg$x.qtde <- length(cfg$x.lim) -1
cfg$y.qtde <- length(cfg$y.lim) -1
cfg$all.qtde <- cfg$x.qtde * cfg$y.qtde
cfg$x <- paste0('X', seq_len(cfg$x.qtde ) )
cfg$y <- paste0('Y', seq_len(cfg$y.qtde ) )

#Sequencias de todos os quadrantes
cfg$all <- paste0('q', seq_len(cfg$all.qtde ) )
names(cfg$all) <- cfg$all

train %>% filter(dropoff_longitude >= cfg$x.min &
                 dropoff_longitude <= cfg$x.max &
                 dropoff_latitude >= cfg$y.min & 
                 dropoff_latitude <= cfg$y.max) -> trainTreino

Monta o mapa, para facilitar vizualização

mapa = matrix( data=cfg$all, nrow = cfg$y.qtde)
colnames(mapa) <- cfg$x.lim[-1]
row.names(mapa) <- cfg$y.lim[-1]

Cria atributos do mapa

mapa.att <- data.frame(names=cfg$all, row.names = cfg$all)
mapa.att <- mutate( mapa.att,
                    x = rep(cfg$x, each = 110),
                    y = rep(cfg$y, 110),
                    x.from = rep( cfg$x.lim[1:110], each = 110),
                    x.to = rep( cfg$x.lim[2:111], each = 110),
                    y.from = rep( cfg$y.lim[1:110], 110),
                    y.to = rep( cfg$y.lim[2:111], 110),
                    x.center = (x.to - x.from) / 2.0 + x.from,
                    y.center = (y.to - y.from) / 2.0 + y.from)

Calculando a qtde. de dropoff no intervalo

QtdePorQuadrante = 0
for (i in 1:12100) {
  trainTreino %>% filter(dropoff_longitude > mapa.att$x.from[i] &
               dropoff_longitude < mapa.att$x.to[i] & 
               dropoff_latitude > mapa.att$y.from[i] &
               dropoff_latitude < mapa.att$y.to[i]) %>% count() -> QtdePorQuadrante[i]
}

#Validação da qtde. de viagens
SomaViagens <- unlist(QtdePorQuadrante)
sum(SomaViagens)
## [1] 266312

Mapa de calor dropoff

cores.mapa  = c("white", "#63be7b", "#ffeb84",  "#f8696b")
gradiente.mapa = colorRampPalette(cores.mapa)(n = 30)
 plot_ly(x = cfg$x,
         y = cfg$y,
         z = matrix(unlist(QtdePorQuadrante), ncol=110, byrow = TRUE),
         colors = gradiente.mapa,
         type = "heatmap")

Plot no Mapa de Nova York

#http://maps.googleapis.com/maps/api/geocode/json?latlng=40.75436,-73.98382&sensor=false

rm(MaiorePontos)
## Warning in rm(MaiorePontos): objeto 'MaiorePontos' não encontrado
MaiorePontos <- filter(mapa.att, (mapa.att$x == 'X27' & mapa.att$y == 'Y23'))
MaiorePontos <- rbind(MaiorePontos, filter(mapa.att, (mapa.att$x == 'X22' & mapa.att$y == 'Y17')))
MaiorePontos <- rbind(MaiorePontos, filter(mapa.att, (mapa.att$x == 'X27' & mapa.att$y == 'Y3')))
MaiorePontos <- rbind(MaiorePontos, filter(mapa.att, (mapa.att$x == 'X60' & mapa.att$y == 'Y26')))

ny.map =get_map(location=c(filter(mapa.att, (mapa.att$x == 'X55' & mapa.att$y == 'Y55'))$x.center,filter(mapa.att, (mapa.att$x == 'X55' & mapa.att$y == 'Y55'))$y.center), zoom = 14,
                source = "google", maptype = "terrain", color = "color")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.75526,-73.98544&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ny.map.2016 <- ggmap(ny.map, base_layer =
                       ggplot(aes(x = x.center, y = y.center), data = MaiorePontos), extent = "device")
## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead
ny.map.2016 + geom_point(size = I(1),colour="red", alpha = 2/3)
## Theme element panel.border missing
## Theme element axis.line.x.bottom missing
## Theme element axis.ticks.x.bottom missing
## Theme element axis.line.x.top missing
## Theme element axis.ticks.x.top missing
## Theme element axis.line.y.left missing
## Theme element axis.ticks.y.left missing
## Theme element axis.line.y.right missing
## Theme element axis.ticks.y.right missing
## Theme element plot.title missing
## Theme element plot.subtitle missing
## Theme element plot.tag missing
## Theme element plot.caption missing

Analises

trainTreino <- mutate(trainTreino, MesPickup = str_sub(trainTreino$pickup_datetime,start = 6,end = 7))

trainTreino <- mutate(trainTreino, SoHoraPickup = str_sub(trainTreino$pickup_datetime,start = 12, end=13))

# Agrupamento por mês, pode-se identificar o mês com maior qtde. de viagens.
AgrupamentoMes <- trainTreino$MesPickup
AgrupamentoMes <- data.frame(AgrupamentoMes)
group_by(AgrupamentoMes) %>% count(AgrupamentoMes) -> Mesprincipal
View(Mesprincipal)

# Agrupamento por hora, para analisarmos a correlação entre a qtde. de viagens por horário, e velocidade média por horário. Conseguimos concluir que a correlação (-0.8893697) é inversamente proporcional, indicando quanto maior a qtde. de corridas no horário, menor é a velocidade média.
trainTreino %>% group_by(SoHoraPickup) %>% summarise(QtdeViagensHora = n(), EVmedia = sum(Vmedia), VMedia = EVmedia / QtdeViagensHora) -> Horaprincipal

View(Horaprincipal)

cor(as.double(Horaprincipal$QtdeViagensHora), Horaprincipal$VMedia)
## [1] -0.8893697
# Agrupamento por carro e por mês, para analisarmos a correlação entre a qtde. de viagens e o valor ganho. Conseguimos concluir que a correlação (0.8958311) é diretamente proporcional, indicando quanto maior a qtde. de viagens, maior o valor ganho.
trainTreino %>% group_by(vendor_id, MesPickup) %>% summarise(QtdeViagensMes = n(), EValor = sum(valor)) -> vendedoresAgrupamento

View(vendedoresAgrupamento)

cor(vendedoresAgrupamento$QtdeViagensMes, vendedoresAgrupamento$EValor)
## [1] 0.8958311
# Verificando também a correlação entre o valor e a duração, concluimos que é diretamente proporcional, e inferimos que o cálculo é baseado basicamente nas durações das viagens.
cor(trainTreino$valor,trainTreino$trip_duration)
## [1] 0.9967498
#cor(AnaliseCorrida2d$valor,AnaliseCorrida2d$distManMilhas)
 
trainTreino %>% filter(trainTreino$valor < 60 & trainTreino$distManMilhas < 2 & trainTreino$SoHoraPickup == 19) %>% select(vendor_id, distManMilhas,trip_duration, valor) -> AnaliseCorrida3d

regressao <- lm(AnaliseCorrida3d$valor~AnaliseCorrida3d$trip_duration, data = AnaliseCorrida3d)
summary(regressao)
## 
## Call:
## lm(formula = AnaliseCorrida3d$valor ~ AnaliseCorrida3d$trip_duration, 
##     data = AnaliseCorrida3d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.5302  -1.0217   0.0839   0.9663  22.7199 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.2847058  0.0421536   77.92   <2e-16 ***
## AnaliseCorrida3d$trip_duration 0.0141592  0.0000448  316.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.668 on 14305 degrees of freedom
## Multiple R-squared:  0.8748, Adjusted R-squared:  0.8747 
## F-statistic: 9.991e+04 on 1 and 14305 DF,  p-value: < 2.2e-16
z <- plot(AnaliseCorrida3d$valor,AnaliseCorrida3d$trip_duration)
grid(z)
abline(regressao)

Modelagem

trainTreino %>% filter(trainTreino$valor < 60 & trainTreino$distManMilhas < 5) %>% select(vendor_id, trip_duration, valor) -> AnaliseCorrida2d

#trainTreino %>% select(vendor_id, valor, distEucMilhas) -> AnaliseCorrida2d

plot(AnaliseCorrida2d[,-1])

rm(modelo.corrida)
## Warning in rm(modelo.corrida): objeto 'modelo.corrida' não encontrado
set.seed(2)
modelo.corrida = kmeans(AnaliseCorrida2d[,-1], centers = 4)
nossoplot <- plot(AnaliseCorrida2d[,-1], col=modelo.corrida$cluster, pch=21, cex=1)

nossoplot + points(modelo.corrida$centers, col = 4:1, bg=1:4, pch = 24, cex =1, lwd = 1)

## integer(0)
trainTreino %>% filter(trainTreino$valor < 60 & trainTreino$distManMilhas < 2 & trainTreino$SoHoraPickup == 19) %>% select(vendor_id, distManMilhas,trip_duration, valor) -> AnaliseCorrida3d

#trainTreino %>% select(vendor_id, trip_duration, valor, distManKM) -> AnaliseCorrida3d

#plot(AnaliseCorrida3d[,-1])

set.seed(2)
modelo.corrida = kmeans(AnaliseCorrida3d[,-1], centers = 4)
plot(AnaliseCorrida3d[,-1], col=modelo.corrida$cluster, pch=21, cex=1)

#points(modelo.corrida$centers, col = 4:1, bg=1:4, pch = 24, cex =1, lwd = 1)
AnaliseCorrida3d %>%
  mutate(cluster = modelo.corrida$cluster) %>% 
  plot_ly(data = .,
        x = ~trip_duration, y = ~valor, z = ~distManMilhas,
        text = ~vendor_id,
        type = 'scatter3d',  mode ='markers',
        color = ~cluster,
        size = rep(1, dim(AnaliseCorrida3d)[1] ), sizes=c(3.0 ) )